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I — i ' 

We propose a new method for the evaluation of the particle density and kinetic pressure profiles in 

#^ ■ inhomogeneous one-dimensional systems of non-interacting fermions, and apply it to harmonically 

confined systems of up to A^ = 1000 fermions. The method invokes a Green's function operator 
in coordinate space, which is handled by techniques originally developed for the calculation of the 
density of single-particle states from Green's functions in the energy domain. In contrast to the 
Thomas-Fermi (local density) approximation, the exact profiles under harmonic confinement show 

>0 , negative local pressure in the tails and a prominent shell structure which may become accessible to 

observation in magnetically trapped gases of fermionic alkali atoms. 

o 

The techniques which have led to the achievement of Bose-Einstein condensation in vapors of bosonic atoms ||1J-U] 
*7^ , are currently being used to trap and cool dilute gases of fermionic alkali atoms 0]. Under magnetic confinement 
t^ ' the s-wave collisions between pairs of fermions in a single hyperfine level are suppressed by the Pauli principle, and 
the geometry of the trap can be adapted to have cylindrical symmetry with a transverse confinement which may be 
hundreds of times stronger than the longitudinal one. It is thus possible to experimentally realize quasi-onedimensional 
C^ , inhomogeneous systems of almost non-interacting fermions at very low temperature and high purity. 
G ' A number of one-dimensional (ID) physical models can be solved exactly H and their solution serves as a test 

I , of approximate theories and contributes to the understanding of real systems. Some important examples are the 
^ ' determination of the ground state and excitation spectrum of a hard-core Bose gas in ID H and the solution of the 
iv Kronig-Penncy model for the electron energy bands in a ID crystal lattice [pi. The spectral and transport properties 

^ of other ID systems of non-interacting electrons have been studied as models for polymers and quantum wires, using 
'~^, Green's function methods [P 12| for which ingenious techniques such as a decimation/renormalization procedure 



^ 



O 
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3|,|14| have been developed. 

^ I In this work we present a new method for the exact evaluation of the ground-state particle density profile in a 

^^^ ■ spin-polarized ID system of up to large numbers of non-interacting fermions in arbitrary spatial confinement. In 
C , essence we show that, just as the single-particle density of states in the energy domain can be obtained by powerful 
Green's function methods, similar techniques yield the particle density in the ID space domain. In fact, our method 
also allows the evaluation of higher moments of the one-body density matrix: we focus here on its second moment, 
p^ which is simply proportional to the kinetic energy density and to the kinetic pressure. 
f^ , As an application of the general method we give results for the particle density and kinetic pressure profiles of a 

^^ ' degenerate Fermi gas in harmonic confinement. This model is directly relevant to the current experiments on atomic 
3 : Fermi gases and we show that the shell structure noticed for the particle density in earlier theoretical studies in 3D 
^ ' |£5|j|6[ is greatly enhanced in ID. We also use our exact results to test the Thomas-Fermi (local density) approximation 
_IL in dependence of the number of fermions in the confined gas. 

a: 
o. 

Q . I. GENERAL FORMULATION 

>■ 

• '^ ■ The one-body Dirac density matrix for a system of N non-interacting fermions at zero temperature can be expanded 

rN ' on the single-particle wavefunctions ^pi{x) = {x jV'i) as p{xi, X2) = J2i=i '4'i{xi)^i{x2)- By using the representation of 
j3 i the translation operator this becomes 

N 

p{x,,X2) = Y.^:{x,)e-^P^'^---^H,{x^) , (1) 

2=1 

showing how distant points are correlated through the momentum operator p. Expansion in powers of the relative 
coordinate r — xi~ X2 yields physical observables such as the particle density profile n{x), 

N 

n{x) = p{x + r/2, x ~ r/2%^^ = ^(^, | 5{x - x,) \ ^,) (2) 



and the kinetic pressure Pix), 






m dr^ ' '' ° 2?7T, 



This is twice the kinetic energy density. 

The main idea of this Letter is to rewrite Eqs. (||) and (||) as the imaginary part of the ground-state average of 
suitable operators related to the Green's function in coordinate space G{x) — {x — x + ie)~^ . We have 



1 ^ 

(x) = hm Im V(V'. \Gix)\ ^,) (4) 



and 



P{x) ^-- hm Im V(V^, \—G{x)\ ^,) . (5) 



TT e^0+ ^ — ' m 

i—1 

G{x) can then be treated by methods analogous to those used for treating Green's functions in the energy domain. 

The equivalence between expressions (^ and 1^ is easily proved in the coordinate representation, where the density 
profile in Eq. (Q) reads 

n{x) = lim Im y^ / dxi\ipi{xi)\'^ — , (6) 

■K e-*o+ ^-^ J X — Xi+ie 

t—i 

yielding Eq. (||) when one takes the limit e ^ 0+. The equivalence between expressions (JSJ) and @ for P{x) is 
similarly proved. 

Evidently, this method can be applied to all ID systems which may be described by single-particle orbitals: one only 
needs to know the representation of the position and momentum operators on such a basis. Hence, interactions could 
also be included in evaluating the particle density through the use of Kohn-Sham single-particle orbitals. Models of 
displacement fields (such as those induced by impurities) may also be studied directly without previous evaluation of 
orbitals. 

II. NON-INTERACTING FERMI GAS IN HARMONIC TRAP 

As already noted, a ID Fermi model is relevant to the spin-polarized fermionic vapors in magnetic confinement [g|, 
where it is possible to realize experimentally a ID configuration by making use of very anisotropic axially symmetric 
traps. At low temperature only the transverse ground state of the trap is populated and the vapor can be described 
by an effective ID harmonic Hamiltonian. 

An analytic expression for the particle density of this system has been given by Husimi |l^,nq] in terms of the 
wavefunction of the 7V-th fermion in the trap. However, a calculation of the density profile and the kinetic pressure by 
his approach is limited to small values of N and has been reported only for TV = 1 and 2 [n8| . Our method allows us 
to efficiently evaluate these ground state properties even for quite large numbers of particles. Of course, the simplicity 
of the representation of the position and momentum operators in this system makes it a favorable example. 

A. Particle density profile 

We start with the calculation of the particle density. For a linear harmonic oscillator the position operator is 
represented as a: = (a + a^)/^/2 on the basis {| ^pi)} of the eigenvectors of the Hamiltonian. As usual, the creation 
and destruction operators satisfy the relations a\tpi) = \/i — 1 1 "04- 1) a-nd a'^lipi) = v'|'0i+i)- A straightforward 
procedure for evaluating the profile (0) is to invert the matrix {x — x + ie) and to calculate its trace on the submatrix 
of its first Nx N block (TrAr). We employ the relation Tr^rQ = d [lndet(Q~^ -I- AIat)] /i9A|a=o, where In is a diagonal 
semi-infinite matrix with its first N eigenvalues equal to 1 and null elsewhere p9| , to obtain 

1 d 

n(x) = lim Im-— - [Indetfx — a; + ie + AIjv)l\_n ■ (7) 

TT e^o+ dX ' 



The calculation of the determinant in Eq. (0) is conveniently performed by the recursive algorithm developed in [gO| 
Renormalization of the tridiagonal operator R = x — XIn allows us to write 



det {x — R + ie) — TT (a; ~ Ofe + is) (8) 

with oi — —A, Ofc+i — —A + ^k(x — Ofe + ie)~^ for 1 < fc < TV and ak+i = ^k{x — at + ie)~^ for k > N. 

The scheme given in Eqs. (7) and (g|) is easily implemented numerically. In practice we have performed the 
calculation of the determinant (8) up to the product of its first M terms, which corresponds to inverting an M x M 
matrix. We have checked the convergence of this approximation by increasing the dimension M and correspondingly 
decreasing the value of e. 

In Fig. |l]we report the density profile n{x) for A^ = 5, 10 and 20 fermions, with M = lO'"^ and e — 0.01. The exact 
profiles are also compared with those given by the Thomas- Fermi approximation (LDA), 

nLDA(a;) = i(2iV-a;2)i/2 (g) 

TT 

(in units such that h = 1, m = 1 and uj = 1). The exact profiles contain N oscillations, which become smaller in 
relative amplitude as N increases. Without any special numerical efforts we have evaluated the exact density profile 
up to iV = 1000: this would otherwise require the calculation of Hermite polynomials up to the 1000*'' degree. For 
large N the oscillations are so small in relative amplitude that their smoothing in the LDA profile becomes reasonably 
accurate, except for the region of the tails. 

B. Kinetic pressure profile 

The particle density profile that we have evaluated above is the analogue of the density of single-particle states in 
the energy domain. In the space domain other single-particle quantities acquire physical interest, as is the case for 
the kinetic pressure in Eq. (ra). Wc show how also for this function one can profitably resort to a renormalization 
technique. 

Taking Q — iPG{x), Eq. (H) can be written as 

1 d 

Fix) ^ lim ln\—-\\ndei{x-x + ie + \lNP^)\. n- (10) 



The matrix [x — x + ie + WnP^) appearing in Eq. (10) is pentadiagonal on the first N rows, owing to the form of the 
operators x and p = i{a^ ~ a)/\/2 in the basis of the energy eigenstates {| f/'i)}- The calculation of the determinant 
of such a matrix can again be performed by the recursive algorithm given in |2Cf| . Renormalization of the operator 
K — X — XInP^ is made on blocks of dimension 2 for the pentadiagonal part and on blocks of dimension 1 for the 
tridiagonal part. This allows us to write 

f n5=r^^' det(x - A, + is) Y\Zn+3{^ - ^k + is) (even N) 
dei{x - K + ie) ^ I (11) 

I n^=i+'^^' clet(x - A, + le) YlZN+ii^ - «fe + *£) (odd N) . 

The renormalized 2x2 blocks Aj satisfy the recursion relation 

Aj ^ Aj + Bjj-i(x - ij-i + ie)-^Bj^i,j (12) 

for j > 1 and Ai = Ai. The matrices Aj, Bjj^i and Bj^ij are submatrices of the operator K, which are defined as 
follows: 



^, -A(2j -3/2)gA,_ 2,+i V(2J - l)/2 , .,„. 

^^ ^ V(2j - l)/2 -A(2j-l/2)0^_2, '' ^ ^ 



Bj,j+i - 



\^23{2j -1)129^-2,+! Q_ \ .^4^ 

^ \^2j{23 + l)/2eN-2j) ^ ^ 



and 

^+''^"1 AV2j(2j + l)/2 0^-2,-2 j' ^ ^ 

with 0fe = 1 for fc > and 9k ^ otherwise. The recursion relation for the elements 5,^ is again ak+i = ^k/ix — ak+ie), 

the first elements being ak+i = ^k{x + ie - [Ak/2h2 - [Ak/2]2i ■ [Ak/2]i2/{x + ie - [Ak/2]ii)}^^ with A; = A^ + 2 for 
even N and fc = A^ + 3 for odd N. We have studied the convergence of the determinant in Eq. (O) as for the case of 
the particle density profile. 

In Fig. the kinetic pressure P{x) is plotted for A^ = 5, 10 and 20 with AI ~ 10^ and e = 0.01, together with the 
profiles PhBAix) evaluated in the Thomas- Fermi approximation, 

Plda = ^i2N ~ x'f/^ (16) 

on 

The exact kinetic pressure shows N oscillations and has the peculiarity of being negative in the tails. This microscopic 
quantum effect, which is missing in the local density description, reflects the fact that in the low density region the 
kinetic energy decreases with increasing density. We have checked that our results agree with those reported in |lj] 
for A^ = 1 and 2, and carried out the calculation of P{x) up to A'' = 1000. The kinetic pressure profile for N=1000, 
as shown in Fig. ra, looks almost indistinguishable from the LDA prediction but stills presents a region of negative 
kinetic pressure in the tails (inset in Fig. 0) . 

In conclusion, in this Letter we have given a general formula for the exact particle density and kinetic pressure 
profiles of a ID many-fermions system in terms of a Green's operator in coordinate space. We have made use of the 
decimation/renormalization procedure and of other recursive techniques, originally developed to evaluate the spectral 
properties of quasi-lD systems in solid state physics, to efficiently calculate the exact density profiles of a harmonically 
confined non-interacting Fermi gas. Within the same general scheme the particle density could also be evaluated by 
employing a suitable Kirkman-Pendry relation [pi| , as will be reported elsewhere. We have verified that for large 
number of atoms (A^ — 1000) the local density approximation reproduces reasonably well the exact profiles except for 
the region of the tails, where the exact kinetic pressure is negative. 

We believe that the present method opens the way for a novel approach to the equilibrium properties of spatially 
inhomogeneous ID systems. The expressions here derived can be extended to finite temperature and to calculate 
partial density profiles for subgroups of atoms. The kinetic energy density functional can be studied through the 
calculation of the function P[x{n)]/2 where x{n) is obtained by local inversion of the exact profile n{x). Pressure 
fiuctuations will become accessible to study through the evaluation of higher moments of the one-body density matrix. 
The density profiles of the harmonically trapped Fermi gas in ID, showing a prominent shell structure as displayed 
in our calculations, could become observable in experiments on alkali vapours. 
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FIG. 1. Exact particle density profile (bold lines) for N = ^, 10 and 20 harmonically confined fermions, compared with the 
corresponding profiles evaluated in the local density approximation. Positions are in units of the characteristic length of the 



harmonic oscillator auo = \/ft/(mtj) and the particle density in units of a 
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FIG. 2. Exact kinetic pressure profile (bold lines) for N = ^, IQ and 20 harmonically confined fermions, compared with the 
profiles evaluated in the local density approximation. Positions are in units of a^o = \/h/{r 
units of hiijar}. 
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FIG. 3. Exact kinetic pressure profile (bold lines) for A'^ = 1000, compared with the profile evaluated in the local density 
approximation. The inset shows an enlarged view of the tail of the profiles. The units are as in Figure 0. In this calculation 
we have employed a matrix of dimension M — 10^ and chosen e — 10~^ (see notations in the text). 



